/*************************************************************************** 

Implimentation of NeedleMan Wunsch  algoritham for 
PairWise Alignment.

Usage : <Needle> <Seq File Name> <Scoring Matrix> <Penalty> <Affinity>
****************************************************************************/
#include <iostream.h>
#include <stdlib.h>
#include <stdio.h>
#include <string>
#include <string.h>
#include <fstream.h>


int readScoringMatrix(char * filename);
int readSequences(char * filename);

char seq[2][1000];	// need to be allocated Dynamically later.

int matrix[23][23];

long **nwMap=NULL;
char **arrowFlag=NULL;

long gap=-8;		// Need to Read From Scoring Matrix 
long affn=-8;

int main(int argc,char* argv[])
{

	if( argc <5 )
	{
		cout<<"Usage : Needle  SeqFile  Matrix  Panalty  Affinity"<<endl;
		exit(1);
	}

	readScoringMatrix(argv[2]);
	readSequences(argv[1]);

	
 ofstream out("out.txt");
	for (int x=0; x<23 ; x++)
	{
		for (int y=0 ; y<23 ; y++)
		{
			out<<matrix[x][y]<<"   ";
		}
		out<<endl;
	}	
out<< "nd of the matrix"<< endl;

	gap=atoi(argv[3]);

	long rowLen=strlen(seq[0]);   //-------------
	long colLen=strlen(seq[1]);   //|
	

	nwMap=new long *[colLen+1];
	for(int i=0; i<colLen+1; i++)
	{
		nwMap[i] = new long[rowLen+1];
	}

	arrowFlag=new char *[colLen+2];
	for( i=0; i<colLen+2; i++)
	{
		arrowFlag[i] = new char[rowLen+2];
	}


//	arrowFalg = new char[colLen*rowLen+2];

	nwMap[0][0]=0;
	 
	for(int j=1; j<colLen+1; j++) // Gap penality in first col.
	{
		//nwMap[j][0] = j*gap;
		  nwMap[j][0] = gap+affn*(j-1);

	}

	
	for(int k=1; k<rowLen+1; k++) // Gap panalty in First row
	{
		//nwMap[0][k] = k*gap;
		  nwMap[0][k] = gap+affn*(k-1);
	}
	

 
 /* for(i=0 ; i<colLen+1 ;i++)  // first col 
  {
	 for(j=0 ; j<rowLen+1 ;j++) // first row
	{
		out<<nwMap[i][j]<<"\t";

	}
		out<<endl;

  }
*/


  for(i=0 ; i<colLen ;i++)  // first col 
  {
    int m=0 , n=0;

	switch(seq[1][i])
	{
	case 'A' :  { m=0 ; break; }
	case 'R' :  { m=1 ; break; }
	case 'N' :  { m=2 ;	break; }
	case 'D' :	{ m=3 ;	break; }
	case 'C' :	{ m=4 ;	break; }
	case 'Q' :	{ m=5 ;	break; }
	case 'E' :	{ m=6 ;	break; }
	case 'G' :	{ m=7 ;	break; }
	case 'H' :	{ m=8 ;	break; }
	case 'I' :	{ m=9 ;	break; }
	case 'L' :	{ m=10;	break; }
	case 'K' :	{ m=11;	break; }
	case 'M' :	{ m=12;	break; }
	case 'F' :	{ m=13;	break; }
	case 'P' :	{ m=14;	break; }
	case 'S' :	{ m=15;	break; }
	case 'T' :	{ m=16;	break; }
	case 'W' :	{ m=17;	break; }
	case 'Y' :	{ m=18;	break; }
	case 'V' : 	{ m=19;	break; }
	case 'B' :	{ m=20;	break; }
	case 'Z' :	{ m=21;	break; }
	case 'X' :	{ m=22;	break; }
	default  :	{ m=23;	break; }
	} // end of the switch
	
  
  for(j=0 ; j<rowLen ;j++) // first row
  {
  	switch( seq[0][j])
	{
	case 'A' :  { n=0 ; break; }
	case 'R' :  { n=1 ; break; }
	case 'N' :  { n=2 ;	break; }
	case 'D' :	{ n=3 ;	break; }
	case 'C' :	{ n=4 ;	break; }
	case 'Q' :	{ n=5 ;	break; }
	case 'E' :	{ n=6 ;	break; }
	case 'G' :	{ n=7 ;	break; }
	case 'H' :	{ n=8 ;	break; }
	case 'I' :	{ n=9 ;	break; }
	case 'L' :	{ n=10;	break; }
	case 'K' :	{ n=11;	break; }
	case 'M' :	{ n=12;	break; }
	case 'F' :	{ n=13;	break; }
	case 'P' :	{ n=14;	break; }
	case 'S' :	{ n=15;	break; }
	case 'T' :	{ n=16;	break; }
	case 'W' :	{ n=17;	break; }
	case 'Y' :	{ n=18;	break; }
	case 'V' : 	{ n=19;	break; }
	case 'B' :	{ n=20;	break; }
	case 'Z' :	{ n=21;	break; }
	case 'X' :	{ n=22;	break; }
	default  :	{ n=23;	break; }
	} // end of the switch
	
// check here for logical errors.
	
	//	nwMap[i+1][j+1]
	//	arrowFalg[1]='D';
		
		long d=matrix[m][n]+nwMap[i][j];

		// to take care of affinity
		

		long s = nwMap[i+1][j]+gap;
		// to take care of affinity

		long u =nwMap[i][j+1]+gap;
		
		char flag;

		if(d>=u)
		{
			flag='D';
			nwMap[i+1][j+1]=d;
		}
		else
		{
			flag ='U';
			nwMap[i+1][j+1]=u;
		}

		if(nwMap[i+1][j+1]>=s)
		{
			;
		}
		else
		{
			flag ='S';
			nwMap[i+1][j+1]=s;
		}

		arrowFlag[i+1][j+1]=flag;
	//	*arrowFalg =flag;
    //        arrowFalg++;

  } // inner for loop.
 
  } // outer for loop.

	//*arrowFalg++ ='\0';
	cout<<"Score :"<<nwMap[i][j]<<endl;


  for(i=0 ; i<colLen+1 ;i++)  // first col 
  {
	 for(j=0 ; j<rowLen+1 ;j++) // first row
	{
		 out<<nwMap[i][j]<<" "<<arrowFlag[i][j]<<"\t";

	}
		out<<endl;

  }

  out<<endl;
  
  
// deleting memory.

	for(i=0; i<colLen+1; i++)
	{
		delete []nwMap[i];
		nwMap[i]=NULL;
	}

	delete []nwMap; nwMap=NULL;

// Finding aligned sequences.

// create two arrays whose length is sum of two sequences.

	char * top= new char[colLen+rowLen+2];
	char * bot=new char[colLen+rowLen+2];
    char *top1=top;
	char *bot1=bot;

//	long rowLen=strlen(seq[0]);   //-------------
//	long colLen=strlen(seq[1]);   //|

j=colLen ;i=rowLen;
while((i>0)&&(j>0))
{

		if(arrowFlag[j][i]=='D')
		{
			*top=seq[0][i-1];
			top++;
			*bot=seq[1][j-1];
			bot++;
			
			i--; j--;

			continue;

		}
		else if(arrowFlag[j][i]=='U')
		{
			*top='-';
			top++;
			*bot=seq[1][j-1];
			bot++;
			
			i; j--;

			continue;
	

		}
		else
		{
			*top=seq[0][i-1];
			top++;
			*bot='-';
			bot++;
			
			i--; j;

			continue;
	

		}
}

/*
if((i==0)&&(j>0))
{
  
  *bot=seq[1][j-1];
	bot++;
	*top='-';
	top++;
	i--; j;
	

}else
if((i>0)&&(j==0))
{

  *top=seq[0][i-1];
	top++;
	*bot='-';
	bot++;
	i--; j;

} */

while(j>0)
{
  
  *bot=seq[1][j-1];
	bot++;
	*top='-';
	top++;
	i; j--;
	

}
while(i>0)
{

  *top=seq[0][i-1];
	top++;
	*bot='-';
	bot++;
	i--; j;

} 



*top='\0';
*bot='\0';

top=strrev(top1);
bot=strrev(bot1);



out<<"\n"<<bot<<endl;
out<<top<<endl;


	for(i=0; i<colLen+1; i++)
	{
		delete []arrowFlag[i];
		arrowFlag[i]=NULL;
	}

// writing the aligned sequence.

char * sym=new char[colLen+rowLen+2];

for(i=0 ; i<strlen(bot); i++)
{

char flag;

if( bot[i]==top[i] )  	                      flag='|';

//X-X

else if( (top[i]=='X') && (bot[i]!='-') )     flag=':';

else if( (bot[i]=='X') && (top[i]!='-') )     flag=':';

//A L I V F M P

else if( ( (top[i]=='A') || (top[i]=='L') || (top[i]=='I') || (top[i]=='V') ||
	       (top[i]=='F') || (top[i]=='M') || (top[i]=='P') )  &&
		   
		 ( (bot[i]=='A') || (bot[i]=='L') || (bot[i]=='I') || (bot[i]=='V') || 
		   (bot[i]=='F') ||( bot[i]=='M') ||( bot[i]=='P') ))

											 flag=':';
// D E K R

else if( ( (top[i]=='D') || (top[i]=='E') || (top[i]=='K') || (top[i]=='R') )&&

		 ( (bot[i]=='D') || (bot[i]=='E') || (bot[i]=='K') || (bot[i]=='R') ))

											 flag =':';
// S T C N Q H Y W

else if( ( (top[i]=='S') || (top[i]=='T') || (top[i]=='C') || (top[i]=='N') || 
		   (top[i]=='Q') || (top[i]=='H') || (top[i]=='Y') || (top[i]=='W'))&&

         ( (bot[i]=='S') || (bot[i]=='T') || (bot[i]=='C') || (bot[i]=='N') || 
		   (bot[i]=='Q') || (bot[i]=='H') || (bot[i]=='Y') || (bot[i]=='W') ))


											 flag=':';

 //B : N D
 
else if( (top[i]=='B') && ( (bot[i]=='N')||(bot[i]=='D')) )   flag=':';

else if( (bot[i]=='B') && ( (top[i]=='N')||(top[i]=='D')))    flag=':';

//Z : E Q 

else if( (top[i]=='Z') && ( (bot[i]=='E') ||(bot[i]=='Q')))   flag=':';


else if( (bot[i]=='Z') && ( (top[i]=='E') ||(top[i]=='Q')))   flag=':';

else flag=' ';


sym[i]=flag;
}
sym[i]='\0';

cout<<sym<<endl;

/*
out<<"#############################################################################"<<endl;
out<<"#"<<endl;
out<<"#	Program :	Needle man wuench"<<endl;
out<<"# Report  :       "<<endl;
out<<"#############################################################################"<<endl;
out<<endl<<endl;
*/
delete []sym;
delete []top;
delete []bot;
	
return 0;
}






// to read blosum matrix
int readScoringMatrix(char * filename)
{
  int i,j;
  char ch , x[25],temp[101]="#";

  ifstream in(filename);
  if(!in)
	  return 1;

  while(temp[0] == '#')	 
	in.getline(temp, 100,  '\n' );

    for(i=0;i<23;i++)
	{
		in>>ch;
		
		x[i]=ch;
		for(j=0;j<23;j++)
			in>>matrix[i][j];
		
		in.getline(temp, 100,  '\n' );

	}

	in>>ch;
	in>>gap;
	cout<<"The Gap dist is "<<gap;
   return 0;
 }

int readSequences(char *filename)
{
//	char seq[2][1000];
	ifstream in(filename);
	if(!in)
	{
		cout<<"No seq file found"<<endl;
		return 1;
	}
	in.getline(&seq[0][0] , 1000, '\n'); // header
	in.getline(&seq[0][0] , 1000, '\n'); // input sequence

	seq[1][0]=' ';
	while(seq[1][0]!='>')
		in.getline(&seq[1][0] , 1000, '\n'); // header

	in.getline(&seq[1][0] , 1000, '\n'); // input sequence

	cout<<seq[0]<<endl;
	cout<<"this is the second sequence"<<endl;
	cout<<seq[1]<<endl;

	return 0;
}






	